Ballistic transport in disordered graphene 
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An analytic theory of electron transport in disordered graphene in a ballistic geometry is devel- 
oped. We consider a sample of a large width W and analyze the evolution of the conductance, 
the shot noise, and the full statistics of the charge transfer with increasing length L, both at the 
Dirac point and at a finite gate voltage. The transfer matrix approach combined with the disorder 
perturbation theory and the renormalization group is used. We also discuss the crossover to the 
diffusive regime and construct a "phase diagram" of various transport regimes in graphene. 



PACS numbers: 73.63.-b, 73.22.-f 
I. INTRODUCTION 

Recent successes in manufacturing of atomically thin 
graphite samples 1 ^' 3 -^ (graphene) have stimulated in- 
tense experimental and theoretical activity 6 -* 7 -. The key 
feature of graphene is the massless Dirac type of low- 
energy electron excitations. This gives rise to a num- 
ber of remarkable physical properties of this system dis- 
tinguishing it from conventional two-dimensional metals. 
One of the most prominent features of graphene is the 
"minimal conductivity" at the neutrality (Dirac) point. 
Specifically, the conductivity 3 ' 4 ' 5 of an undoped sample 
is close to e 2 /h per spin per valley, remaining almost con- 
stant in a very broad temperature range — from room 
temperature down to 30mK. 

Several recent theoretical works addressed transport 
in disordered graphene samples. It was found that lo- 
calization properties depend strongly on the nature of 
disorder— iVi 11 ' 12 ! 13 ! 14 which determines the symmetry 
and topology of the corresponding field theory. The lo- 
calization is absent provided a certain symmetry of clean 
graphene Hamiltonian is preserved in the disordered sam- 
ple, see Ref. [T3- 

One possibility is that disorder preserves a chiral sym- 
metry of massless Dirac fermions. This situation is, 
in particular, realized when the dominant disorder is 
due to corrugations of graphene sheet (ripples) and/or 
dislocations.-^ 5 - The conductivity in such chiral-symmetric 
models has been shown 1 ^ to be exactly e 2 /nh (per spin 
per valley) at the Dirac point. While being temperature 
independent, this value is, however, less by a factor of 
~ 3 than the experimentally measured values. 

Another possibility, a long-range randomness, was 
studied in Ref. [ID. This type of disorder does not 
mix the two valleys of the graphene spectrum, which 
leads to emergence of a topological term in the corre- 
sponding field theory (unitary or symplectic a- model). 
The peculiar topological properties protect the sys- 
tem from localizatio n 11 ! 12 ! 13 ! 16 ! 17 . It is worth mention- 
ing that a topologically protected metallic state emerg- 



ing in graphene with long-range random potential also 
arises at a surface of a three-dimensional Z-i topological 
insulator. 18,19 

A number of numerical simulations of electron trans- 
port in disordered graphen o 17 ' 20 ! 21 ' 22 ! 23 ' 24 confirmed the 
absence of localization in the presence of long-range ran- 
dom potential. The main quantity studied numerically 
in most of these works is the conductance G of a finite- 
size graphene sample with a width W much larger than 
the length L. This setup allows one to define the "con- 
ductivity" a = GL/W even for ballistic samples with L 
much shorter than the mean free path I. Remarkably, in 
graphene at the Dirac point, such ballistic "conductivity" 
has a universal value e 2 /irh in the clean case i 25 ' 26 This 
setup was studied experimentally in Refs.|27. 28. 29,3(3 and 
the ballistic value e 2 /nh was indeed observed for large as- 
pect ratios. This geometry of samples is particularly ad- 
vantageous for the analysis of evolution from the ballistic 
to diffusive transport. 

A complete description of the electron transport 
through a finite system involves not only the conduc- 
tance but also higher cumulants of the distribution of 
transferred charge. The second moment is related to the 
current noise in the system. The intensity of the shot 
noise is characterized by the Fano factor F. For clean 
graphene, this quantity was studied in Ref. [2(| Surpris- 
ingly, in a short and wide sample (W ^> L) the Fano fac- 
tor takes the universal value F = 1/3, that coincides with 
the well-known result for a diffusive metallic wirei^I This 
is at odds with usual clean metallic systems, where the 
shot noise is absent (F = 0). The Fano factor F = 1/3 
in clean graphene is attributed^ 6 , to the fact that the cur- 
rent is mediated by evanescent rather than propagating 
modes. Furthermore, the whole distribution of transmis- 
sion eigenvalues for the massless Dirac equation in a clean 
sample with W ^> L at the Dirac point agrees with that 
of mesoscopic metallic wires in the diffusive regime^ 

The effect of disorder on the shot noise was studied 
numerically in Refs. |23||24| . where the value of the Fano 
factor F w 0.3 was found across the whole crossover 
form ballistics to diffusion. The Fano factor close to 1/3 



was also observed at the Dirac point experimentally] 29 i 30 
When the chemical potential was shifted away from the 
Dirac point, the Fano factor decreased, then showed 
an intermediate shoulder at F w 0.15, and finally ap- 
proached zero for largest gate voltages (carrier concen- 
trations). 

While both diffusive and clean limits have been ad- 
dressed analytically, only numerical and experimental re- 
sults for the intermediate regime of ballistic transport 
through disordered samples have been available so far^ 
The aim of this paper is to fill this gap. We develop 
the analytic theory of electron transport in disordered 
graphene in the ballistic geometry (L <C W, I) and cal- 
culate the full statistics of the charge transfer for both 
zero (the Dirac point) and large concentration of car- 
riers. We also discuss the crossover to diffusive regime 
and construct the overall "phase diagram" of transport 
regimes. 

The structure of the paper is as follows. We begin in 
Sec. [II] with the introduction of the model and derivation 
of a general transfer-matrix equation. In Sec. IIIII we cal- 
culate transport properties of a clean sample. In Sec. IIVI 
the disorder is included in the lowest order of the pertur- 
bation theory. The resummation of leading higher-order 
corrections to the counting statistics is performed within 
the renormalization group approach in Sec. [V] For the 
case of random potential we present an evidence in fa- 
vor of a universal scaling of the distribution function of 
transmission coefficients valid at the Dirac point for sam- 
ples of arbitrary size (covering both ballistic and diffusive 
regimes). We summarize the results and discuss the per- 
spectives in Sec. I VII Technical details are relegated to 
Appendices [5J and [C] 



II. TRANSFER-MATRIX TECHNIQUE 

We start with introducing our model and the general 
formalism of transfer matrix technique. For graph ene, 
this approach was employed in Refs. 23.25.26.34.35.36. 

We will adopt the single- valley model of graphene. 
More specifically, we will consider scattering of electrons 
only within a single valley and neglect intervalley scat- 
tering events. Indeed, a number of experimental results 
show that the dominant disorder in graphene scatters 
electrons within the same valley. First, this disorder 
model is supported by the odd-integer quantization^ 4 ^ 
of the Hall conductivity, a xy = (2n + l)2e 2 /h, represent- 
ing a direct evidence 1 ^ in favor of smooth disorder which 
does not mix the valleys. The analysis of weak local- 
ization also corroborates the dominance of intra-valley 
scattering 3 ^. Furthermore, the observation of the linear 
density dependence^ of graphene conductivity away from 
the Dirac point implies that the relevant disorder is due 
to charged impurities and/or ripples . 10 i 20 i 38 i 39 i 40 i 41 Due 
to the long-range character of these types of disorder, the 
intervalley scattering amplitudes are strongly suppressed 
and will be neglected in our treatment. Finally, apparent 
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FIG. 1: (Color online) Schematic setup for two-terminal 
transport measurements. Graphene sample of dimensions 
L x W is placed between two parallel contacts. We assume 
W ^> L throughout the paper. 

absence of localization at the Dirac point down to very 
low temperatures 3 -" 4 ^ can be explained only by some spe- 
cial symmetry of disorder. The most realistic candidate 
model is the long-range randomness which does not scat- 
ter between valleys^ 

The single-valley massless Dirac Hamiltonian of elec- 
trons in graphene has the form (see, e.g., Ref. 0) 

H = v ap + V(x,y), V(x,y) = a^{x,y). (1) 

Here <7 M (with /i = 0, x, y, z) are Pauli matrices acting on 
the electron pseudospin degree of freedom corresponding 
to the sublattice structure of the honeycomb lattice, <x = 
{cr x ,<jy\, and the Fermi velocity is vq ~ 10 8 cm/s. The 
random part V(x, y) is in general a 2 x 2 matrix in the 
sublattice space. Below we set H = 1 and vo = 1 for 
convenience. 

We will calculate transport properties of a rectangu- 
lar graphene sample with the dimensions L x W. The 
contacts are attached to the two sides of the width W 
separated by the distance L. We fix the a; axis in the 
direction of current, Fig. [TJ with the contacts placed at 
x = and x = L. We assume W ^S> L, which allows us to 
neglect the boundary effects related to the edges of the 
sample that arc parallel to the x axis (at y — ±W/2). 

Following Ref. [26, metallic contacts are modelled as 
highly doped graphene regions described by the same 
Hamiltonian ((T|). In other words, we assume that the 
chemical potential Ep in the contacts is shifted far from 
the Dirac point. In particular, Ep 3> e, where e is the 
chemical potential inside the graphene sample counted 
from the Dirac point. (All our results are independent 
of the sign of energy, thus we assume e > through- 
out the paper.) A large number of propagating modes 
exists in the leads, all belonging to the circular Fermi 
surface of radius pf = Ep/vq. These modes are la- 
belled by the momentum p n = 2nn/W in y direction 
with \n\ < N — Wpf/2it. Particular boundary condi- 
tions at y = ±W/2 shift the quantized values of p n by 
a constant of order 1/W. However, this constant has no 
significance in the limit W 3> L when many channels 
participate in electron transport. 

Clearly, the transverse momentum p n is preserved in 
the clean system. We will use the mixed momentum- 
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coordinate representation, with the wave function 4'„(x) 
bearing a vector index n in the space of transverse mo- 
menta supplemented by a 2-spinor structure in pseu- 
dospin (sublattice) space. The eigenstates of the clean 
Hamiltonian Hq = VQtrp have the direction of pseudospin 
parallel to the electron momentum. It is convenient to 
perform the unitary rotation^ in the pseudospin space 
if) = £^ with C = (cr x + a z )/y/2 which transforms a x 
to the diagonal form: La x C) = a z . Hence the two com- 
ponents of the rotated spinor correspond to right- and 
left-propagating waves^ 2 - ip = {ip R , tp L }- In terms of the 
new function ip n (x), the Schrodinger equation H'fy = evp 
acquires the for m 23 i 34 



dx 



(<7xPn + icr z e)tp n 



E 



(2) 



The matrix U nm {x) represents the operator CV(x,y)C^ 
in the mixed momentum-coordinate representation 



(3) 



A standard description of electron propagation in- 
volves the scattering matrix S. This is a unitary matrix 
relating the amplitudes of incident and outgoing waves 



^(0) 
^ R (L) 



S 



S 



r t' 
t r' 



(4) 



The elements t, t' and r, r' are matrices in channel space 
formed by transmission and reflection amplitudes, re- 
spectively. The unitarity condition S^S = 1 ensures con- 
servation of particle number. 

A closely related formulation is based on the transfer 
matrix T which expresses the waves at the point x = L 
through those at x = 0: 



ip R (L) 



= T 



il> R (0) 

V> L (o) 



T = 



r't' 



-t'-\ t'- 1 



(5) 

This description is convenient due to the simple mul- 
tiplicativity property: T{x%, X2)T(x2, x±) = T{x^,x\). 
The current conservation is provided by the identity 
T^a z T = a z . 

By definition, the transfer matrix T(x2 1 x\) yields a 
solution to the Schrodinger equation ^ in the form 
VK 3 ^) = T(x2, Xi)tp(xi). Transfer matrix itself, as a 
function of its first argument, obeys the same Schrodinger 
equation with the initial condition T(x, x) = 1. In a 
clean sample the solution depends only on the difference 
x 2 — x\ and is diagonal in channel space: 

T^l(x 2 ,xi) = S nm e^q>[(a x p n + ia z e)(x 2 - xi)]. (6) 

In order to include disorder as a perturbation, it is 
convenient to cast the Schrodinger equation @ into an 
integral form. In terms of transfer matrix the integral 



equation reads 

T(X2,X 1 )=T^{X 2 ,X 1 ) 

,x)a z U(x)T(x,xi). (7) 



The transport statistics of the sample is expressed in 
terms of transmission eigenvalues T n — the eigenvalues 
of the matrix tH. One can extract these transmission 
eigenvalues from the upper left element of the transfer 
matrix ([S]). The first two moments of the transferred 
charge distribution determine the conductance (by Lan- 
dauer formula) and the Fano factor— 



4e 2 

G=— Tr(t% 



F = 1 — 



Tr(tH) 2 
Tr(itt) 



(8) 



The factor 4 in the expression for the conductance ac- 
counts for the spin and valley degeneracy. 



III. CLEAN GRAPHENE 

We will first analyze transport properties of a clean 
graphene strip. In the "short and wide" geometry (W ^> 
L) we are considering, the total number of channels par- 
ticipating in charge transfer is large. This allows us to 
replace summation over channels by integration. From 
now on, we will identify channels by the dimensionless 
momentum p = p n L in y direction and integrate over 
this momentum according to 



Y^wjdp 

^ T. I 2tt 



(9) 



The transfer matrix T^°\ and hence its upper- left 
block , are diagonal in channels. Using the explicit 
form of the clean graphene transfer matrix, Eq. ((6|), one 
calculates the transmission eigenvalues 3 ^ 



1 



p 2 sinh 2 y/p 2 — (eL) 



P 



{eLf 



(10) 



For the conductance and Fano factor we obtain from 



Eq. 



G = 



2e 2 W 
irhL 



dp T p , 



F=l- 



JdpT } 



J dp T p 



(11) 



The result of numerical integration of Eq. (jlip is shown in 
Fig. [21 A detailed analytical analysis of the two limiting 
cases of small and large energies is presented below. 



A. Transmission distribution and counting 
statistics 

It is convenient to introduce the distribution function 
P(T) of transmission eigenvalues (TTUj) . This distribu- 
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FIG. 2: Energy dependence of the (a) conductance and the 
(b) Fano factor of the clean sample with W S> L. Solid lines 
show numerical results. Low energy asymptotics Eq. (|20l) is 
plotted by dashed lines while dotted lines correspond to high 
energy limit Eqs. (|33[1 and (|34p . Asymptotical curves provide 
a very good approximation to the exact result in the whole 
range of energies. 



tion function provides a measure in the space of chan- 
nels which is, by definition, equivalent to the integration 
measure ©, 



P(T) dT 



Wd\p\ 



(12) 



According to (fTUl) , there is one-to-one correspondence be- 
tween the transmission eigenvalue < T < 1 and the 
absolute value |p| of the momentum; an extra factor of 2 
in the right-hand side of Eq. (fT2"|) accounts for the double 
degeneracy between channels with momenta p and —p. 

Very generally, the distribution function P(T) deter- 
mines the full statistics of the charge transfer. Specif- 
ically, one defines the counting statistics k(x) = 
J2 N e ixN V(N), where V(N) is the probability that N 
particles are transferred within a measurement time in- 
terval t m . Then \nn(x) is the generating function for 
cumulants, 



mK (%) = 



k\ 



((N k )). 



(13) 



It can be related to P(T) in the following way 4 ^ (we 
assume zero temperature and retain the factor of 4 taking 
into account the spin and valley degeneracy in graphene) : 



4e 2 f 
In «( X ) = —Vt m J dTP(T) ln(l - T 



+ e lx T), (14) 



where V is the applied voltage. In particular, the first 
two of the cumulants ((N k )) determine the conductance 
G and the shot noise power S via ((N)) = (N) = Vt m G 
and ((N 2 )) = Vt m S. According to Eqs. Q3J), ([13]). one 
has 

4e 2 f 4e 2 f 

G=— dTTP(T), S=— dTT(l-T)P(T), 

(15) 

and the Fano factor 



F= £ = 1 JdTT 2 P(T) 
G f dTTP(T) ' 



(16) 



The relations (HU), (fT5")) . (HHJ) are of general validity and 
equally applicable to the clean and disordered system. 
All the information about scattering, both at the inter- 
face with leads and in the bulk of the system, is encoded 
in the transmission distribution P(T). Clearly, Eqs. (|15p . 
(I16p for the conductance and the shot noise are equivalent 
to Eq. ®. 



B. 



Low energies: ei < 1 



In the low energy limit, we calculate the distribution 
function P(T) in the form of a power series in the small 
parameter eL. In order to perform this calculation, we 
first invert the function T p given by Eq. (|10[) keeping 
terms of the second order in eL: 



p(T)=p (T) 



(eL) 



1 



Po(T) pl{T) 



Po (T) — arccosh —== . 



Now we substitute this expression into Eq. 
tain the distribution 



(17) 
(18) 

W and ob- 



P(T) = 
IT' 



W dp(T) 
~^L 

1 



dT 



2-kL T^/T^T 



1 + {eLf 



1 + T 



VT^T 

PUT) 2pl(T) 



(19) 



It is worth noticing that by definition J dTP(T) should 
give the total number of open channels Wpf/tt in the 
leads. In fact, the logarithmic divergence at T — * of the 
normalization of Eq. (fIT)|) is cut off at the lowest trans- 
mission eigenvalue T m j n ~ exp(— 2ppL). This small-T 
cutoff is, however, immaterial for the calculation of the 
moments (conductance, noise, etc). 

At zero energy, the function P(T) reproduces the well- 
known Dorokhov result 43 for a diffusive wire. This is, 
in particular, the reason for the 1/3 Fano factor in 
graphene. 2 - The fact that the clean graphene sample is 
characterized by exactly the same form of the transmis- 
sion distribution as a generic diffusive wire is highly non- 
trivial. We will show below (Sees. IIVI and |V| that this 



remarkable correspondence remains valid in the ballistic 
regime when leading disorder effects are incorporated. 

Using the distribution (|19|) . we obtain the following 
results for the conductance and the Fano factor of clean 
graphene at low energies, eL < 1, 



G 



Cl 



f'2 



0.101, 



35C(3) 124C(5) 

3tt 2 V a 

28C(3) 434C(5) 4572C(7) 

157T 2 7T 4 7T 6 



c 2 (eL) 2 ], 

(20) 

(21) 

-0.052. (22) 



At the Dirac point (e = 01 . Eq. (|20|) reproduces the earlier 
analytical results of Refs. |25||261 Low energy asymptotics 
is shown with dashed lines in Fig. [5] 




FIG. 3: Integration contour used in Eq. (|29|l . 



dpf(p) 



cL 



>, du 



x [f(eWl~u2) + f(-eWl-u 2 )]. 



(26) 

Transforming Eq. ([24]) to the new variable ([25]) and 
averaging over oscillations (see Appendix [5] for details), 
we obtain 



C. High energies: eL > 1 

When the Fermi energy e in the sample is far from 
the Dirac point, many conducting (T ~ 1) channels are 
opened. In this regime, the conductivity and higher mo- 
ments of the transmission distribution are essentially lin- 
ear in e with small oscillating corrections (see Fig. [2]). 
These oscillations are due to interference effects: conduc- 
tance is relatively enhanced and the noise is suppressed 
when a channel exhibits resonant transmission with T 
close to 1. This phenomenon is similar to the Fabry- 
Perot resonances. 

We begin with the calculation of the main (propor- 
tional to e) part of the transmission distribution func- 
tion and will return to the oscillatory correction later in 
this section. It is convenient to find first the generating 
function of transmission moments, defined as 



T(z) = J2 « n-1 Tr(tU) n = Tr 



(23) 



T{z) 



We 



ttVI 



ir 



du 



irV-l — zu^ 
We K(z) 

IT 



E(z) 



z 



VT 



(27) 



Here K{m) and E{m) are complete elliptic integrals of 
the first and second kind^ with the parameter to. 

The function T(z) is regular at the point z = 0. The 
coefficients of the series expansion near this point pro- 
vide the moments of transferred charge distribution, see 
Eq. (|2"3")) . The transmission distribution function P(T) is 
related to T(z) by the linear integral equation 



P(T) 



dT = T(z), 



(28) 



which follows from Eq. (|23|) . In order to solve this equa- 
tion for P(T), we note that the function T{z) has a 
branch cut along real axis running from 1 to oo. We in- 
tegrate Eq. I|28p along a contour going from z = 1/T—iO 
to z = 1/T + iO encircling the point z = 1, see Fig. [3] 
This integration yields 



This function appears to be very useful for the forthcom- 
ing calculation of the transport properties of a disordered 
sample. In this section, we apply it to the clean system. 
According to Eqs. (f2"3"| and (fTU|) , we have 



T( \ W [ d P 



l-z 



1 ■ 2 

p sin 



-i -i 



{eLf - p 2 



(24) 

The integrand oscillates rapidly in the interval — eh < 
p < eL. This interval of momenta contains all open 
(well-conducting) channels and thus provides the main 
contribution to the generating function. At high ener- 
gies, it is convenient to introduce a new variable u, such 
that 



V =Pr, 



eLy/T 



(25) 



1/T+iO 



dz 
2ttI 



1/T+iO 



Hz) 



P(f)df 



dz 



1/T-iO 



1/T-iO 



2m(T- 1 



P(f)df. (29) 



To find the distribution function, we calculate the deriva- 
tive of the above equation with respect to T and obtain 



P(T) = 



jT + iO) - T(l jT - % 0) 
2niT 2 ' 



(30) 



This identity establishes a relation between the distri- 
bution function P(T) and the jump of T{z) across the 
branch cut at the point z = l/T. In other words, equa- 
tion (|28|) solves the corresponding Riemann-Hilbert prob- 
lem. 
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To find the explicit formula for the distribution P(T), 
we perform an analytic continuation of the expression 
([27| from the vicinity of the point z = to z = 1/T±i0 
and substitute the result into Eq. (|30[) . This yields 



We K(T) - g(T) 

1 J ^ 2 ryr^r ' 



(31) 



This distribution function provides a full transport de- 
scription at high energies (up to small oscillatory cor- 
rections discussed below). We note that Eq. (|3Tj) does 
not take into account almost closed (evanescent) chan- 
nels with p n > e (and thus T -C 1). Estimating their 
contribution, we find that it is suppressed by a factor 
(eL)~ 4 compared to the main term, thus yielding a neg- 
ligible contribution to the charge transfer. 

Using the above distribution (or, equivalently, the gen- 
erating function), we calculate the asymptotics of the 
conductance and Fano factor in the high-energy limit 
eh > 1, 



G 



We, 



F 



(32) 



We recall that when passing from Eq. ([24f to Eq. ([27f , 
we neglected the oscillatory contributions to the gener- 
ating function. A more accurate calculation accounting 
for these oscillations is presented in Appendix lAl [see Eq. 
(|A6[) ]. The results for the conductance and the Fano fac- 
tor read 



G 



F = - 
8 



We 



1 - 



sin(2eL-7r/4) 
+ 20F(eL) 3 /2 
9sin(2eL - tt/4) 



20F(eL) 3 / 2 



(33) 
(34) 



These results are in a good agreement with the high- 
energy behavior of G and F calculated numerically, see 
Fig.H 

Let us emphasize that transport properties of the sys- 
tem at high energies depend on the particular model of 
the contacts i 46 ! 47 In our calculation we assume that the 
boundaries between graphene and the leads are sharp. 
This model is well justified if the actual extension d of 
the transitional region at the interface is small compared 
to the electron wavelength in graphene. This wavelength 
is energy dependent: it tends to infinity at e — and 
decreases with increasing e. Thus the small-energy re- 
sults, Sec. IIIIB) are universal and not influenced by the 
microscopic details of the interface provided the size of 
the boundary transitional region is much smaller than 
the length of the sample L. On the other hand, the re- 
sults of the current section are applicable for not too 
high energies, e <C 1/d. For higher energy, the electron 
wavelength becomes comparable to d and the transmis- 
sion properties of the sample become nonuniversal. In 
the extreme high-energy limit, ed 3> 1 the boundary be- 
comes adiabatically smooth. This, in particular, leads 
to the vanishing Fano factor because the semiclassically 



propagating electrons are either transmitted or reflected 
without any uncertainty. 

Our results for the energy dependence of the conduc- 
tance and the Fano factor in clean graphene are in agree- 
ment with the findings of Refs. l26ll34l . where the sum over 
transmission channels was evaluated numerically for a fi- 
nite (but sufficiently large) ratio W/L. Experimentally, 
such a ballistic setup was studied in Refs. 1291130 . Most 
of the experimental observations reasonably agree with 
our results. The "conductivity" GL/W (which is equal 
to 4e 2 /-7r/i at the neutrality point, as expected for a bal- 
listic sample) increases roughly linearly with energy e. 
The Fano factor has a value close to 1/3 at the Dirac 
point and decreases when one moves away from the Dirac 
point, showing a tendency to saturate at F « 0.15, which 
is not far from the value 1/8 we have obtained in the 
high-energy regime. Measurements on other samples re- 
veal that very far from the Dirac point the Fano factor 
decreases again, reaching a value as low as 0.02. Appar- 
ently, the intermediate plateau corresponds to the high- 
energy regime L~ l e d^ 1 investigated in our work, 
while the vanishing of the Fano factor at still higher elec- 
tron concentrations corresponds to the ultra-high-energy 
range, e 3> d -1 . It is not quite clear to us why the oscilla- 
tory structures are not observed in experimental data. A 
possible explanation is that the length L of the sample in 
the experiment varies as a function of the y coordinate, 
leading to a suppression of the oscillations. 



IV. INCLUDING DISORDER: PERTURB ATIVE 
TREATMENT 

So far, we have considered the transport properties of a 
clean graphene sample. In the present section we include 
disorder on the level of the leading perturbative correc- 
tion. As discussed in the beginning of Sec. [ill we neglect 
the intervalley scattering. Further, we will assume the 
Gaussian statistics for disorder components U M that de- 
termine the random part V(x,y) of the Hamiltonian ([1]) 
acting within a single valley, 

(W(x,y))=0, (35) 
(V»(x, y)V(x', y')) = -x',y- y'). (36) 

This type of randomness is realized when the scattering 
is due to impurities in the substrate separated by a thick 
(compared to the lattice constant) clean spacer layer from 
the graphene plane. The intervalley matrix elements of 
the disorder potential are then exponentially suppressed 
and can be safely neglected. A more realistic case of long- 
range charged impurities with 1 jr potentials can also be 
treated perturbatively within the Gaussian model, but 
with an energy-dependent scattering amplitude.— We 
will briefly discuss modifications of the results in the case 
of Coulomb- type impurities in Sec. IV Fl 

The correlation function w^{x,y) is even with respect 
to both arguments and is peaked at short (compared to 
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the wavelength in the sample) distances, being hence al- 
most a delta function. At the same time, we will have 
to keep a small but non-zero correlation length in or- 
der to regularize ultraviolet singularities arising at the 
intermediate stage of our calculation. The results of the 
calculation will not depend on this correlation length. 

In the transfer-matrix approach, it is convenient to 
convert the correlation function to momentum represen- 
tation in y direction. The x dependence of can be 
safely replaced by the delta function without generating 
any singularities. Thus we introduce new dimensionlcss 
functions a M according to 



Sjxl 
W 



(37) 



The functions a^(q) vary slowly with q. They are al- 
most constant at low values of momentum and decay at 
the large scale of inverse correlation length. We will ex- 
press the transport characteristics of the system in terms 
of four constants 



MO). 



(38) 



These parameters are nothing but the amplitudes of the 
effective delta functions in Eq. ([36]) . 



(39) 



and correspond to the intravalley scattering parameters 
used in Ref. [n] (there it was assumed that a x = a y = 
a x /2). 

In the present section we will calculate the first disor- 
der correction to the transport properties of a graphene 
sample. Specifically, we will find a linear-in-a^ contribu- 
tion to the function P(T). It is convenient to introduce 
short-hand notations for inverse transmission amplitudes 
and probabilities of the clean sample, 



K = 1/4°), H n = \h n \ 2 = i/r(°). 



(40) 



Further, it will be useful to label the types of disorder 
(p = 0, x. y, z) by a pair of binary indices £, r) = ±, ac- 
cording to 



OLy = a , 



(41) 



We develop the perturbative expansion by iteratively 
solving Eq. (JT]). Then we single out the upper- left block 

of the matrix T(L, 0) thus obtaining P . Up to the 
second order in V the result is 



(f^ )mn $mnh n + A mn S mn h n S A T 



(42) 



Here A m „ is the linear correction to the transfer matrix, 



A = -i 



dx 



T^(L,x)a z U(x)T^(x,0)} , (43) 



where the subscript 1, 1 refers to the upper- left block in 
the right/left-mover space, see Eq. ([5]). 

The last term in Eq. (|42p represents the contribution of 
the second order in disorder amplitudes U{x). Since we 
are interested in the correction to transport coefficients of 
the linear order in (and thus quadratic in U), we can 
perform disorder averaging of this term using Eqs. (|36p . 
(I37p . Then the integration over ^-coordinate in this term 
is trivial due to the delta function in the correlator (137)) 
and the multiplicativity property of the transfer matrix. 
We have also used the relation dx5(x) = 1/2. [This 
identity holds because the delta function in Eq. (|3"T|) is 
a replacement for some symmetric sharply peaked func- 
tion.] As a result, we get 



A n 



nL 
W 



^2 £ a £,ri(lr,iL - q n L). 



(44) 



The sum over intermediate states I in the last term of 
Eq. (|42[) converges due to a non-zero correlation length 
of disorder, encoded in the momentum dependence of 
in Eq. gH) . 

Now we substitute the expression (|42l) and its Hermi- 
tian conjugate into Eq. (|23[) and then expand T(z) up to 
the second order in A and first order in A. Performing 
the disorder averaging of terms containing A, we obtain 
the following expression for the generating function T{z): 



n 



i 



+ (H n + z)C„ 



{H n -zf (H n -z) 2 (H m -z) 



with 



(45) 



(46) 



In a general case, the two matrices B mn and C mn are 
very complicated functions of m and n. We will simplify 
further analysis by considering two limiting cases of low 
and high energy. 



A. Low energies: eL <JC 1 

We have already calculated the lowest order correction 
to the distribution function P(T) due to small energy [see 
Eq. (fT9| ], Now we are going to find the lowest disorder 
correction at exactly zero energy. To the main order, 
these two contributions merely add up. 

At zero energy, there are no propagating modes in 
graphene, all the channels are evanescent. In this sit- 
uation, it is convenient to use the transverse momentum 
p = p n L instead of index n to label the channels accord- 
ing to Eq. ((9l). The bare transfer matrix {6} at e = 
simplifies to 



q-(o) _ f coshp sinhjA 
p I sinhp coshp J 



(47) 
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The quantities h n and H n introduced in Eq. (|40|) take 
the form 

h p = cosh p, H p = cosh 2 p. (48) 
The matrix A nm defined by Eq. ([43]) now becomes 

A pq = J dx$y-iVp q {x) cosh[p — (p + q)x/L] 

— iV p x q (x) cosh[p — (p — q)x/L] 
+ V p y q (x) wah\p-(p-q)x/L] 

+ iV; q {x) sinh[p - (p + «)as/L]}. (49) 

Two types of averages [Eq. (|46|) ] arise in the calculation 
of transport properties. These averages are the result of 
applying Eq. (f3"6")) to the product of two A pq matrices 
and subsequent integration over single [due to the delta 
function in Eq. (|3T|) ] position x. This yields 



Br, 



nL 
W 



h p h q ^2 a ^,{p-q) 

t,v=± 



w 



£ cosh(p — rjq) 
sinh(p + -qq) 

P + m y 

sinh 2p + rj sinh 2q 



(50) 



2(p + m) 



(51) 



Now we substitute Eqs. (pHj). fgQ]l. and (JSJ) into Eq. 
5]) and separate the resulting expression into four parts, 



(52) 
(53) 



T{z) = T {z) + T^z) + T 2 {z) + T 3 {z), 
W f dp 

W f 
= ~2~kL ^ I d P d 1 ^ a Zn(P ~ l) 

HpH q tanhp tanh q 



zf{H q - z) 



T 2 {z) 



W 



(54) 



(55) 



£,17 

sinh 2p 



dpdq 
P + 7?g 

77 sinh 2q 



{Hp-zf 



(56) 



The first part, J-'q, originating from the first term of Eq. 
(|45[) . is the generating function for the clean sample, 



W arcsiny^ 
W = -L^fz=W 



(57) 



It corresponds to the distribution function P(T), Eq. 
OH), at e = 0. 



The other three terms are disorder-induced corrections. 
The integral in Eq. (|54[) would not be absolutely conver- 
gent if we replace by constants. For this reason, we 
have to retain the momentum dependence of a^ v (orig- 
inating from finite correlation length of disorder) in the 
integrand. Performing first the integration over p + q and 
then over p — q, we get 



£.17 



(58) 



Note that the value of T\ does not actually depend on 
the precise form of the functions a^lp — q), but only on 
their values at p — q = 0. Indeed, the integral over p + q 
and the subsequent integral over p — q are convergent 
even with constant a^ v . The finite disorder correlation 
length is needed only to ensure the absolute convergence 
of the q integral in Eq. (|54[) . 

The integral in Ti, Eq. (|55|) . is absolutely convergent in 
both variables. This allows us to neglect the momentum 
dependence of a^ v (p — q) and replace it by a constant 
from the very beginning, yielding 



T 2 {z) = T${z)Y^i°t-tir)- 



(59) 



The last term is also absolutely convergent in both 
variables. When writing Eq. (|56p . we have simplified the 
integrand by means of symmetrization with respect to 
p <-> r\q. The integrand in Eq. (|56[) can be rewritten in 
the form of a total derivative, 



d 



d 



dp ^ dq 



Hp — H q 



(p + m )(Hp-z)(H q -z) : 



and hence 



T 3 {z) = 0. 



(60) 



(61) 



Collecting all the terms, we finally obtain the following 
generating function 



T(z) = (1 + 2a 



W arcsinTi 
^ ] -L^z=¥- (62) 



We see that the random vector potential, a x<y , does not 
influence transport characteristics of the system in the 
lowest order. In fact, any vector potential is unable to 
alter conductance or higher moments of charge transmis- 
sion at zero energy. We will give a general proof of this 
statement in Sec. IV CI below. Another manifestation of 
this property was found in Refs. 10.49.50] where it was 
shown that the random vector potential does not change 
the conductivity of an infinite graphene sample at the 
Dirac point. 

The distribution of transmission eigenvalues follows 
from Eq. ^1 with the help of identity flH|). Together 
with the energy correction from Eq. (|19p . P(T) acquires 
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the form 



into Eq. (|45|) and averaging over oscillations, we find 



P(T) 



W 



1 



2ttL Ty/T^T 



1 + 2(a - a z 



l + T 



(arccoslf'l I A T) 2anr<>slr(l/\ T) 



(63) 



Remarkably, the functional dependence P(T) is not 
changed by disorder at e = 0. We will discuss the conse- 
quences of this fact in Sec. IV Bl 



We 



du 



7T Jo V /(1- M 2)( 1 _ zu 2) 



6 l^^io V(l-« 2 )(l-^7 



(68) 



The first term in the square brackets gives the generating 
function of the clean sample [Eq. (f2"T|) ]. while the second 
term represents the leading disorder-induced correction, 
^dis(^)- Evaluating integrals in Eq. (|68|) . we express this 
correction in terms of elliptic integrals: 



B. High energies: ei > 1 

The transport properties of a clean graphene sample 
at high energies were considered in Sec. IIII CI The main 
contribution to the conductance and to higher moments 
is proportional to eL and comes from the band of fully 
opened channels with \p n \ < e. In the present section 
we will calculate the disorder-induced correction coming 
from the same channels. As we will show below, the 
relative correction is of the order of a^eL. Since all the 
momentum integrals will be restricted to \p n | < e, we do 
not need the ultraviolet regularization and can neglect 
the momentum dispersion of from the very beginning. 

As appropriate for high energies (see Sec. IIII C|) . we 
will label the channels by variables u and v related to 
p n and p m according to Eq. (125| . In this representation 
the quantities h n and H n introduced in Eq. (|40|) take the 
form 



117., 

•w^) = - 2(1 _ z) E £ - E ^ 

x[(l-Zz)K(z)-E{z)]. (69) 

Expanding this generating function at z — 0, we read- 
ily calculate disorder corrections to the conductance and 
Fano factor. Combining these corrections with the re- 
sults for clean sample, Eqs. (f3"3"| and (f3~4"|) , we obtain 



h u — cos(weL) — i 
H u = cos 2 (ueL) + 



sin(weL) 

U 

sin 2 (ueL) 



(64) 
(65) 



The matrix A uv [Eq. (|4"3")) ] and the averages B uv and 
C uv [Eqs. (|50p and l[5ip] contain rapidly oscillating terms. 
The integration over u and v will average out these oscil- 
lations. For this reason, we can drop all the terms in B uv 
and C uv that are proportional to odd powers of sin(uei) 
or sin(weL), already before calculating the integrals in 
Eq. (|45|) . Furthermore, we discard the contributions that 
are odd functions of p n and/or p m , which corresponds to 
dropping odd powers of vl — u 2 and vl — v 2 . The ma- 
trices B uv and C uv simplify to 



B„ 



ttL 
~W 



E 



£,H U H V 



1 



(66) 
(67) 



G 



We 



sin(2eL-7r/4) 
2^F(eL) 3 / 2 

7T 

— — eL(a + ol x + 3a y + 3a z ) 



(70) 



F = 



9sin(2eL - tt/4) 
2v^(eL) 3 / 2 

+ — eL(3a n 



3a x + l'3a y + 13a z ) 



(71) 



We see that at high energies any disorder suppresses con- 
ductance and enhances noise at the level of the lowest 
perturbative correction. 

To find the disorder correction to transmission distri- 
bution function (|A9|) we perform the analytic continua- 
tion of .Fdis(z) from the vicinity of the point z = to 
z = l/T±iO and apply Eq. The result is 

WLe 2 ^ 
P dis(T) = __2^a Ci; 



2tt 2 



AE(T)[K{\ - T) - £g(l - T)\ + 7r(C - 1) 
1-T 



(72) 



There is, however, a subtlety in determination of Pdi S (T), 
which is related to the singularity of J^disi 2 ) &t z — 1, see 
discussion of a similar problem in the clean case (Ap- 
pendix [XJ. As a result, the distribution function (|72|) 
cannot be applied in the vicinity of T = 1. Specifically, 
we have to impose the bound 



Substituting these expressions together with Eq. (T4~4"l) 



1 - T » {a^eLf. 



(73) 
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At 1 — T ~ (a^eL) 2 , disorder-induced correction ([72]) 
becomes comparable to the main (clean) term [Eq. PTj) ] 
and our perturbative expansion breaks down. It should 
be stressed that this peculiarity in the behavior of P(T) 
near T = 1 does not affect the evaluation of the moments 
using the generating function T{z), which is based on the 
behavior of the latter in the vicinity of z — 0. Indeed, the 
disorder-induced correction Eq. (|69p to T(z) [and thus 
to the moments, Eqs. (|70| and (fTTj) ] is controlled by the 
small parameter a^eL <C 1. Note that at sufficiently high 
energies eL > l/a M , disorder correction to P(T) becomes 
comparable to the clean result in the whole range of T. 
This implies a crossover to the diffusive regime, where the 
perturbative approach developed in the present section 
fails, see Sec. |V1 



V. RENORMALIZATION GROUP AND 
OVERALL PHASE DIAGRAM 

In the previous section we have calculated the lowest 
disorder correction to transport properties of a ballistic 
graphene sample. In the present section we will discuss 
the resummation of higher-order contributions. 

The second-order and all higher terms contain loga- 
rithmic divergences and thus become important when 
system is still in the ballistic regime, L <C I- 
These logarithms are intrinsic for two-dimensional Dirac 
fermions subjected to disorder and were extensively 
studied in various contexts using renormalization group 
tcchniquej 48 i 49 ' 51 i 52 ' 53 i 54 Application of such a renormal- 
ization group (RG) to disordered graphene was developed 
in Refs. THE 

The RG deals with the two-dimensional action describ- 
ing disordered Dirac fermions, 



d 2 x 



>a Al (V>cr M V) 2 



(74) 



where tp and tp are two-component fermionic (anti- 
commuting) fields. The field-theoretical description of 
a disordered system involves also some tool to get rid 
of diagrams with closed fermionic loops. This is usually 
either supersymmetry or replica trick. In both cases, 
the fields acquire additional structure in supersymmet- 
ric or replica space. Equivalently, on can derive the RG 
equations by simply discarding all diagrams that contain 
fermionic loops, without extending the fields. 

The one-loop diagrams contributing to renormalization 
of energy and disorder couplings are shown in Fig. 2J The 
solid lines are propagators of free electrons, 



G<°>(p) = 



<xp 



(75) 



Dashed lines stand for disorder correlators 2tt ^2 a^a^ ® 
a^. We cut the logarithmic divergence in the one- loop di- 
agrams by the running scale parameter A, which has the 



(a) 



(b) 



(c) 



v 



\ / 

A 
/ \ 

(d) 



dimension of length, and obtain the beta function; 



FIG. 4: One-loop diagrams for (a) electron self energy and (b- 
d) scattering amplitude. These diagrams yield RG equations 
ESI) - 173. 



,10.55 

(76) 
(77) 
(78) 
(79) 



dao 

d In A 

da x 



2(a + a z )(a + a x + a y ), 
da,, 



d\nA <91nA 

da z 



— 2ana 



0«z, 



SlnA 

8e 
SlnA 



2(a + a z )(-a z + a x + a y ), 
e(a + a x + a y + a z ). 



Bare values of energy and disorder couplings, which 
are the initial conditions for RG equations, correspond 
to the scale of the order of lattice spacing or disorder 
correlation length. This scale plays the role of ultraviolet 
cutoff in our theory. We will denote it a. After renor- 
malization procedure we obtain renormalized values of 
the parameters at the scale A and also a new effective 
bandwidth 1/A. 

The renormalization proceeds until one of the following 
events happens: (i) the running scale A reaches the sys- 
tem size L, (ii) one of the disorder couplings becomes of 
the order unity, or (iii) the renormalized energy reaches 
the bandwidth. We will discuss these three possibilities 
for particular disorder types below. Once the renormal- 
ization has been performed, we can calculate observables 
by simply applying the perturbation theory. The results 
of previous section for transport characteristics thus re- 
main applicable with bare parameters replaced by their 
renormalized values. 



A. Random scalar potential 

We start the discussion of various disorder types with 
the case of random scalar potential. Let us first con- 
sider the zero energy limit when the only parameter of 
the model is the disorder coupling a$. In the single- 
parameter case, the RG beta function is universal, i.e. 
does not depend on the regularization scheme, within 
the two-loop accuracy. A discussion of the universality 
and the derivation of the second-loop contribution is pre- 
sented in Appendix [Bj The two- loop RG equation reads 



da o 2, 3 
cUnA =2a ° + 2a °- 



(80) 



The disorder strength, quantified by cto, increases in 
course of renormalization. The renormalization process 
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should be stopped when the renormalized value of ao 
becomes of order of unity, so that the perturbative ex- 
pansion of the beta function fails. The corresponding 
scale is the zero-energy mean free path, which we denote 
Iq. To find this length, we express A as function of ao in 
Eq. (|80|) and integrate from the initial value of ao to 1. 
This yields 



l = a^/a^e 



l/2a 



(81) 



The universality of the two-loop equation is evident from 
Eq. (|81[) . The first loop contribution determines the ex- 
ponential factor in Iq while the second loop gives ^fa® in 
the pre-exponent. This parametric dependence of Iq can 
not depend on the regularization scheme. On the other 
hand, the third loop would fix the numerical prefactor in 
Eq. (|8ip. However, the value of the ultraviolet length a is 
itself defined only up to a number within the framework 
of the linearized Dirac Hamiltonian model. 

At scales shorter than the mean free path lo, the renor- 
malized value of ao is given by 



a (A) = 



1 



21n(7 /A) +lnln(/ /A)' 



(82) 



As long as L <C lo [and thus a${L) <C 1], we can describe 
the transport properties by the distribution function (|63|) 
with the renormalized value olq{L) and e — 0. It is worth 
noting that the lowest order perturbation theory used 
for derivation of Eq. (|63|) in combination with the RG re- 
sult (|82|) provides the best possible accuracy within the 
framework of disordered Dirac Hamiltonian. Specifically, 
the second-order terms in the perturbative expansion of 
P(T) in powers of ao(£) would generate the contribution 
of the same order as that of the third-loop correction to 
Eq. (|5D|) . The latter, however, depends on the regular- 
ization scheme and hence is nonuniversal, as discussed 
above. 

A small but non-zero energy does not change the qual- 
itative behavior of the system, as long as the RG flow is 
terminated by the system size. We refer to this situation 
as "ultraballistic regime" . The energy gets renormalized 
according to Eq. (|79p , which is universal only in the one- 
loop order (see Appendix |SJ ■ Using the result (|52")l , we 
solve the RG equation for energy and obtain 



s(A) = 



v /a [21n(/o/A) + lnln(/ /A)] 



(83) 



It is worth mentioning that the renormalized coupling 
and the renormalized energy (|83p are related via 



e 2 (A) a (A) 



a 



(84) 



The value e(L) is to be substituted into Eq. (|63|) along 
with the renormalized value of ceo(L). This yields the full 
description of transport properties for the system in the 
ultraballistic regime. In particular, the conductance and 



the Fano factor are 
1 + 



irh L 



2a +c 1 {eLf 



? = \ 



1 



ao[21n(Zo/i) + lnln(/o/L)] 



a [21n(Z /L) + lnliuV£)] 



(85) 
(86) 



with the constants given by Eqs. (|2"Tj) and (f2"2"|) . 

When the initial (bare) value of energy is increased, 
the renormalized energy eventually becomes comparable 
to the effective bandwidth 1/A before the running scale 
A reaches L (and still before the disorder coupling ao(A) 
reaches unity). The length scale at which e(A) = 1/A 
plays the role of the effective Fermi wavelength A. (In- 
deed, in the absence of disorder, energy is not renormal- 
ized and A = 1/e.) Using Eq. (f8"3"|) , we find 



(87) 



A = -\/2a ln(e/7), 



where 7 is the characteristic disorder-induced energy 
scale 
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a /l = Ae 



-l/2a 



A = 1/a is the initial bandwidth of the model, and we as- 
sumed that e ^> 7. For e < 7, the role of the wavelength 
is played by the mean free path lo- Note that Eqs. ([57]) 
and ([55]) have the same two-loop accuracy as Eqs. (fSTj) 
(|83|) : in particular, the absence of a double logarithm 
term in Eq. (|87p and of ao in the pre-exponent in Eq. 
(155)) is fully controllable. According to Eqs. fST]). (154"]) . 
the renormalized values of the coupling constant and the 
energy at the scale of the wave length are given by 



e(A) 



^2a ln(e/ 7 ) 



a (A) 



21n(e/ 7 )' 



(89) 



In Fig. [5] we show the phase diagram of various trans- 
port regimes. If e < 7 and L <§C Iq or, alternatively, e > 7 
and L <C A, the renormalization terminates by the sys- 
tem size, A = L, and the system is in the ultraballistic 
regime discussed above [see Eqs. (|85| and (|86| ], If e ^> 7 
and A<L<I, the renormalization stops at A = A and 
the running scale does not reach L. We refer to this case 
as "ballistic regime" , since the system size is still smaller 
than the mean free path Z, 



I = 



X 



7ra (A) 



[21n(e/ 7 )] 3/2 . 



(90) 



This valued of the mean free path corresponds to the 
imaginary part of the electron self-energy calculated in 
the Born approximation with renormalized coupling con- 
stant ao(A). Note that for the model with random scalar 
potential, the transport mean free path, which deter- 
mines diffusion coefficient, is twice longer, i tr = 21. 

In the ballistic regime, the renormalized energy is such 
that e(A)L — L/X ^> 1. This means that we have to use 
the high-energy results of Sec. HVBl In particular, with 
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diffusive (antilocalization) 




ultraballistic 



0.4 r 



0.3 



0.2 



0.1 



0.0 L 



FIG. 5: Schematic "phase diagram" of various transport 
regimes in the graphene sample with random scalar potential. 
The lines indicate crossovers between corresponding regimes. 
The shortest sample exhibits ultraballistic transport with the 
conductance and Fano factor given by Eqs. 1)850 and (|86|) re- 
spectively. When the length of the sample exceeds Fermi wave 
length (|87[). ballistic results (|91|l and (|92[) apply. In a sample 
longer than the mean free path (|90[) . diffusive regime estab- 
lishes with the Drude conductivity (|93[) an d the Dorokhov 
distribution of transmission eigenvalues (|94|l . The conductiv- 
ity experiences symplectic antilocalization in this case. 



the renormalized parameters, the conductance and the 
Fano factor, Eqs. (|70|) and (|7Tj) . become 



1 



sin(2L/A-7r/4) 
_ 1+ 2V^(i/A) 3 / 2 
9sin(2L/A-7r/4) 



L 



3L 
4/" 



(91) 
(92) 



In the expressions (|91j) and (|92|) there are two corrections 
to the leading term. The first (oscillating) correction 
exists in the clean limit and is small provided L » A. 
The second correction due to disorder is small only if 
L -C I. This imposes the natural upper bound on the 
ballistic regime: if the system size exceeds the mean free 
path, electron transport becomes diffusive. In this case, 
the system is naturally characterized by the conductivity 
<7, which determines the conductance via the Ohm's law, 
G = aW/L. The Drude expression for the conductivity 
reads 8 



4e 2 8e 2 

= ln( e / 7) . 
7r/iao(A) irn 



(93) 



The distribution function of transmission eigenvalues in 
the diffusive regime is the same as in a usual quasi-one- 
dimensional metallic sample^ 3 - 



fJ 



W 



2-kL tVT^t' 



(94) 



with the dimensionless conductivity g — (7r/i/4e 2 )<7. Tak- 
ing into account interference effects leads to L depen- 
dence of g in this formula, as we are going to discuss. 



FIG. 6: Unified scaling function (|96|l for both ultraballistic 
and diffusive regimes at zero energy in the case of the random 
potential disorder. 



B. Single parameter scaling for random potential 
at zero energy 

Remarkably, the transmission distribution function at 
zero energy appears to be the same in ultraballistic and 
diffusive limits. In both cases it has the form of Dorokhov 
distribution (|94"|) with the parameter 



l + 2cto(L), ultraballistic, 

— k cr(L), diffusive, 
. 4e^ 



(95) 



which has the meaning of the dimensionless conductivity. 
In the ultraballistic regime, the scaling of g is induced via 
renormalization of ao according to equation (|80[) while in 
the diffusive limit, g 3> 1 acquires antilocalization correc- 
tions characteristic for a disordered system of symplectic 
symmetry. This allows us to infer a unified scaling law 
covering both limiting cases 



3 In g 
dlnL 



(5~1) 2 - 
1 



(9 -If 



g - 1 « 1, 
9 > 1. 



(96) 



The scaling function (f9"6"|) is depicted in Fig. [51 It is qual- 
itatively similar to the numerical results of Ref. [22]. 

The applicability of the Dorokhov distribution to the 
diffusive system in the considered geometry requires a 
comment. The original derivation given by Dorokhov 43 
assumes a diffusive system of a quasi- ID geometry (a 
thick wire), with W -C L, On the other hand, our ge- 
ometry is entirely different, W 3> L. This difference is, 
however, of minor importance for the statistics of charge 
transfer as long as the system is a good metal. Indeed, 
there exist alternative derivations of the Dorokhov statis- 
tics that are based on the semiclassical Green function 
formalis m 56 ! 57 , on the sigma-model approach 5 - or on the 
kinetic theory of fluctuations 5 ^ and do not require any 
assumption concerning the aspect ratio of the sample. 
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We stress that the ultraballistic asymptotics of the 
beta function in Eq. is only valid for Gaussian white- 
noise statistics of random potential. The interpolation 
between the two asymptotics of the beta function in Fig. 
[6] implicitly assumes a smooth crossover between the two 
regimes (in particular, without any intermediate fixed 
points), as suggested by numerical simulations ! 17 ' 22 ' 23 

The scaling function characterizes the evolution 
of the dimensionless conductivity with increasing L. 
Whether the full distribution function P(T) retains ex- 
actly its form (parameterized by g only) in the 
crossover remains an open question. Strictly speaking, 
what we know at the moment is that this form of P(T) 
emerges (i) in the clean limit, (ii) in the ultraballistic 
regime within the first order in ao(L), and (iii) in the 
diffusive regime. On this basis, one could speculate that 
this might be an exact statement for the whole crossover. 
Clearly, this is only a hypothesis that requires further 
verification. 



C. Random vector potential 

Let us now consider the situation when the only disor- 
der in the sample is the random vector potential (char- 
acterized by the couplings a x and a y ). This situation is 
physically realized when disorder is due to random cor- 
rugations of the graphene sheet (ripples). The one- loop 
RG equations for random vector potential read^ 9 . 



da,, 



da,. 



= 0, 



ainA din A 



(97) 
(98) 



In fact, the beta function for the disorder couplings is 
identically zero in all loops ; 49 ' 54 i.e., the random vector 
potential is not renormalized. Since the couplings do not 
change with growing system size, the energy follows a 
power-law: 



e(A) 



'A 
a 



(99) 



where a±_ := a x + a y . For not too high energies, the 
RG flow terminates by the system size L (ultraballistic 
regime), so that 



e(L)L = eL(L/a) ai - < 1. 



(100) 

As demonstrated in Sec. lIVl at zero energy the lowest- 
order perturbative correction to the transport coefficients 
is absent in the case when the only disorder is vector po- 
tential. Now we present a general argument showing that 
any given configuration of the vector potential A(x, y) 
does not affect transport properties of the system at zero 
energy. 

The zero-energy Dirac equation takes the following 
form in the presence of vector potential, 



We fix the gauge by requiring that VA = in the bulk of 
the sample and normal component of A vanishes at the 
boundary of the sample. This gauge is widely used in 
the theory of superconductivity and is referred to as the 
London gauge in that context. In this particular gauge 
we can express vector potential using a scalar function 
(j>{x,y) as 



Ax — Tx ) 

ay 



Ay ~ dx 



(102) 



The boundary conditions allow us to fix = at the 
edges of the sample. The function <j> is related to the 
magnetic field B = d x A y — d y A x by the Poisson equa- 
tion V 2 (/> = —B. The existence of a solution to such an 
equation follows from an equivalent electrostatics prob- 
lem: finding the potential of the charge distribution with 
a given density inside a metallic cavity. 

Now we do a pseudo-gauge transformation introducing 
the new wave function according to \P = e a '^^. For 
this new function, the Dirac equation becomes 



T *Vp* = 0, 



(103) 



which is equivalent to the free Dirac equation with no 
magnetic field. The boundary conditions of the London 
gauge fix <j> = outside of the sample. Thus W = # in 
the leads. The transfer matrix of the whole system, and 
hence all the transport properties, is not influenced by 
the vector potential. This result was obtained in Ref. [60 
for a particular configuration of vector potential. Re- 
cently, we have become aware of an alternative proof 
of the general statement by Titov.— The immunity of 
the transport properties to the vector potential holds de- 
spite the fact that the random vector potential problem 
represents a critical theory^ 9 - with the multifractal wave 
function ^(x,y) and a spectrum of multifractal expo- 
nents governed by the disorder strength a± (for review 
see Ref^ 2 -). 

We will use the results of Sec. lIIII obtained for the clean 
sample. Disorder, however, affects the transport proper- 
ties through the disorder-dependent renormalization of 
energy. Thus, in the ultraballistic regime, we can use 
equations (fl"9|) and (|20|) of Sec. IIII 51 with eL replaced by 
its renormalized value given by Eq. (|100|) . 

When the system size L is larger than the Fermi wave- 
length A, the renormalization of energy stops by the 
bandwidth. The value of A is found from the equation 
1/A = e(A), yielding 



A 



Hi) 



a±/(l+a±) 



(104) 



er(p- A)* = 0. 



(101) 



The i-independent renormalized energy e(A) is thus 
given by e(A) = e {A/e) a±/{1+a±) . 

To calculate the transport coefficients in the ballistic 
regime, we substitute the product e(X)L = i/A » 1 
along with the couplings a x and a y into Eqs. (|70|) and 
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(f71~j) . which yields 



r _ e 2 W 
II 



sin(2£/A-7r/4) ' 
1+ 2^(W /2 -4A^ +3a «) 



F = - 



9sm(2i/A-V4) ttL 
1 " 2^F(L/A)3/ 2 + 4A (3 ^ + 13a »> 



Using Eq. (|104p we find the mean free path 



1(e) = 



A 



7raj_e 



£ \ ai/{l+«i) 

A/ 



(105) 

) • 

(106) 

(107) 



The system becomes diffusive when L > I. In contrast to 
the case of scalar potential discussed in Sec. IV Al for ran- 
dom vector potential there is no direct crossover between 
the ultraballistic and diffusive regimes (at zero energy 
the mean free path diverges). 

At finite energy the system belongs to the unitary 
symmetry class. The corresponding field-theory pos- 
sesses a topological termjii which drives the system to 
the quantum-Hall critical point with the universal value 
of conductivity a = (0.5 -j- 0.6)4e 2 //i. Since the dis- 
order coupling aj_ stays non-renormalized, the Drude 
conductivity^ is given by the Born approximation with 
the bare value of a± , 



2e 2 



nhaj_ 



(108) 



and does not depend on energy. The criticality is 
achieved only at very large scales L ~ £ cor . The 
quantum-Hall correlation length £ cor is of the order of 
the unitary-class (second-loop) localization length £ oc 
exp(g 2 ) at which all states would be localized in the ab- 
sence of the topological term, 



1(e) e 



l/4«i 



(109) 



The phase diagram for the random vector potential (Fig. 
[7]) contains four regions: ultraballistic (0 < L < A), bal- 
listic (A < L < I), diffusive (I < L < £ cor ), and critical 

(L > £cor)- 



D. Random mass 

Let us now discuss the transport properties in a sit- 
uation when disorder is modeled solely by the random 
mass term (a z coupling). We are not aware of physical 
realizations of such disorder in graphene. Nevertheless, 
we will consider this case for the sake of completeness. 

The RG equations for the random mass read 



da z 
din A 



= -2ai 



de 
9 In A 



(110) 



It is worth mentioning that there is an exact (valid in 
all loops) relation^ between the beta function (3 Z for the 



critical 

diffusive (weak localization) 



ballistic 
ultraballistic 



L = A 



A 



FIG. 7: Schematic "phase diagram" for the sample with ran- 
dom vector potential. Ultraballistic and ballistic regimes are 
similar to the case of random scalar potential. However, the 
time-reversal symmetry is broken and the system exhibits 
weak (second-loop) localization in the diffusive regime, which 
eventually drives it to the critical state characteristic for the 
quantum Hall transition. There is no disorder-induced energy 
scale in the case of random vector potential, thus we plot the 
phase diagram for all energies up to the band width A, which 
is exponentially larger than 7 in Fig. [5] In the case of random 
mass, the phase diagram is qualitatively the same. 



random-mass coupling a z and (3q for the random poten- 
tial atj: 



Pz(a z ) = ~(3o(-a z ). 



(Ill) 



However, we do not need the two-loop result for the ran- 
dom mass, since a z (A) decreases with growing A and 
therefore the second-loop term never becomes important. 
The one-loop RG equations for a z and e are solved by 



a, (A) 



1 + 2a z ln(A/a) ' 



5(A) = 



a z / „ , A 

— 7TT = (\ 1 + 2a 2 ln-. 

a z (k) V a 



(112) 



(113) 



Thus, while a z (A) decreases with growing length scale 
A, the energy e(A) becomes larger. Energy reaches the 
bandwidth at the Fermi wave length scale A determined 
by e(A) = 1/A. This yields 



A 



1 



1 



e y/l + 2a z ln(A/e)' 



e(A) = e^/l + 2a z ln(A/e). 
At this scale the disorder coupling becomes 

az[ >~ l + 2a z lnA/e 



(114) 
(115) 

(116) 



To describe the transport properties in the ultraballis- 
tic regime (L <C A), we substitute Eqs. (|112[) and (|113|) 
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taken at A = L into Eq. (|63[) and find 
4e 2 W f 2a z 



G = 



Trh L 
■ Cl {eLf[l 



1 - 



l + 2a z \n(L/a) 
2a z ln(L/a)]}, 



{l + c 2 (eL) 2 [l + 2a z ln(L/o)]}, 



(117) 
(118) 



When the system size exceeds the wavelength, L ^> A, 
we use the values for a z and e given by Eqs. (|116[) and 
(j!15[) . The system becomes diffusive when 



L>1 = 



A 



1 



1 



7ra z (A) na z e 
In the ballistic regime A <C I < ( we get 



A 

2a z In — 

e 



e 2 W 



1 - 



sin(2L/A - tt/4) 3L 
1 + 20F(i/A) 3 / 2 ~ IT 
9sin(2L/A-7r/4) 13L 
20F(i/A) 3 / 2 + "IT 



(119) 

(120) 
(121) 



where A is given by Eq. (|114p . 

At zero energy, the system belongs to the supercon- 
ducting symmetry class D (see, e.g., Ref.l52p. Finite en- 
ergy drives the system to the unitary symmetry class A, 
similarly to the case of random vector potential. Again, 
weak (second-loop) localization leads to the quantum- 
Hall critical point at very large scales L > £ cor . The 
Drude conductivity^ is given by the Born approxima- 
tion with the renormalized value of a z from Eq. (|116p . 



4e 2 
3nh 



2 In 



(122) 



The corresponding quantum-Hall correlation length is 
then given by 



1(e) exp 



1/1 „, A 
h 2 In — 

9 \a z e 



(123) 



Thus the case of random mass disorder is very similar 
to the case of random vector potential. We have four 
regimes: ultraballistic (0 < L < A), ballistic (A < L < I), 
diffusive (I < L < £ cor )) and critical (L > £ C or)- The 
schematic phase diagram is the same as for the random 
vector potential case presented in Fig. [7J There is no 
direct crossover between the ultraballistic and diffusive 
regimes. This is related to the fact that at zero energy 
the system is ultraballistic for arbitrary length L. 



E. Generic disorder 

Finally, let us consider the case of generic disorder, 
when all disorder couplings are present. In fact, even if 
only two of the three coupling constants cyq, a±, and a z 



are present at the initial ultraviolet scale a, the third one 
always becomes non-zero with growing system size^ see 
Eqs. flTJJ), dZZ]), and ([78]). The system belongs to the 
unitary symmetry class at all energies and falls into the 
quantum Hall universality classj 11 ' 49 Physically, this sit- 
uation is realized in graphcne when, e.g., both long-range 
vector (ripples) and scalar (charged impurities) potential 
are present j 11 ^ 2 ' 13 

As discussed in Appendix [Bl when two or more 
coupling constants are nonzero, the second-loop beta- 
function becomes non-universal. Therefore, we will deal 
here with the one-loop RG equations. The solution of 
the set of coupled RG equations ([75]). (|77j) . and (J7HJ) is 
analyzed in Appendix [C] It turns out that when the 
initial values of the couplings are of the same order, af- 
ter renormalization the coupling ao (corresponding to 
the scalar potential) dominates. In particular, at zero 
energy the renormalization stops at the scale Iq when 
tto^o) ~ 1' while l ne other two couplings are still much 
smaller than unity (suppressed by a logarithmic factor), 
a±(lo),a z (l n ) ~ (9/8)| lnao| _1 , where ao 1 is the bare 
value of the couplings. 

The phase diagram for the case of generic disorder (Fig. 
[5]) contains four regimes, similarly to the cases of random 
vector potential and random mass: in addition to the 
ultraballistic, ballistic, and diffusive regimes, there is a 
regime of quantum-Hall criticality. On the other hand, at 
variance with the random vector potential and random 
mass problems, the diffusive regime consists of two sub- 
regimes (with weak antilocalization and weak localization 
corrections to the Drude conductivity, respectively). In- 
deed, as discussed above, at the border of the diffusive 
regime (L ~ I) the dominant coupling is ao, which cor- 
responds to the symplectic symmetry class £2- At larger 
scales (L > l c ~ io| lnaol -1 / 2 for e = 0), the gap in the 
Cooperon modes due to the couplings a± and a z becomes 
important^ and only diffuson modes remain, restoring 
the unitary symmetry and leading to the second-loop lo- 
calizing correction to the conductivity. When the renor- 
malized conductivity drops down to the value of the order 
e 2 /h, the critical quantum Hall regime sets in. It is also 
worth mentioning that, similarly to the case of random 
scalar potential, at lowest energies (including e = 0) the 
ultraballistic regime crosses over directly into the diffu- 
sive/critical regime. 



F. Additional comments 

Throughout this paper, we have considered a some- 
what idealized theoretical model. Specifically, we have 
neglected (i) intervalley scattering, (ii) momentum de- 
pendence of scattering amplitude characteristic for scat- 
terers with 1 jr potentials (like charged impurities or rip- 
ples) , and (iii) electron-electron interaction. We are now 
going to discuss, on the qualitative level, a possible in- 
fluence of these effects on our results. 
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FIG. 8: Schematic "phase diagram" for the case when more 
than one disorder type is present in the system. Random 
scalar potential ao becomes dominant in the course of ballistic 
renormalization (see Appendix [C} ; therefore, the ultraballis- 
tic, ballistic, and lowest part of the diffusive regime are similar 
to the diagram Fig. [5] Once the diffusion is established, the 
antilocalization starts but it proceeds only till the length l c at 
which the time-reversal symmetry breaks down.— At longer 
scales the system falls into the unitary symmetry class and ex- 
hibits weak (second-loop) localization. At exponentially long 
scale £ cor the quantum Hall critical state establishes. 

(i) Intervalley scattering. As discussed in the begin- 
ning of the paper, the fact that the dominant disorder 
scattering in experimentally studied graphene samples is 
of intra-valley nature follows directly from the observed 
anomalous, odd-integer quantum Hall effect. The domi- 
nance of the intra-valley scattering also explains why the 
localization is not observed at the Dirac point down to 
very low temperatures. Therefore, the model of decou- 
pled valleys considered in this work is not only of theoret- 
ical interest but is also directly relevant to experiments. 

Still, in any realistic system some amount of inter- 
valley scattering will be present, so that it is natural to 
ask what its influence will be. The weakness of the inter- 
valley scattering implies that the corresponding mean 
free path k n t OT is much larger than the mean free path 
I (induced by the intra-valley scattering and considered 
in the paper). This means that the scale Winter is gener- 
ically located far in the diffusive (or critical) regime. 
The results in the ultraballistic and ballistic regimes, as 
well as in a parametrically broad window in the diffu- 
sive and critical regimes, remain essentially unaffected 
by the inter-valley scattering. At very large distances, 
L 'inter) the intravalley scattering will strongly affect 
the behavior, generically inducing the localization (ex- 
cept for a special case of chiral disorder at the Dirac 
pointiS). 

(ii) l/r impurities. Most of realistic candidates 
for long-range scatterers in graphene samples, such as 
charged impurities and ripples, are characterized by l/r 
potentials. As was shown in Ref. [l(| there is no ballistic 
RG for this type of scatterers; the scale- and energy- 
dependences of disorder-induced effects in the ballistic 



regime are governed simply by the energy dependence of 
the cross-section of an individual scatterer. With these 
modifications, all the considerations in our paper remain 
applicable. In particular, all the phase diagram remain 
qualitative unchanged; one should just use the appropri- 
ate values of the wave length A and of the mean free path 
I. 

(Hi) Electron- electron interaction. The effect of el- 
ectron-electron interaction on the system of disordered 
Dirac fermions constitutes in general a very complex 
problem. In the clean case, the interaction induces 
a logarithmic correction to the velocity^ that can be 
treated within an RG scheme similar to the ballistic dis- 
order RG used in this work. In the disordered case, 
a unified ballistic RG emerge o 65 i 66 describing renormal- 
ization of disorder couplings and of the interaction. In 
Ref. [H corresponding one-loop RG equations are derived 
for time-reversal-invariant disorder and in the limit of 
large number of valleys (simplifying the theoretical treat- 
ment). One can use this interaction-modified RG values 
for renormalizcd couplings entering our results in the ul- 
traballistic and ballistic regimes; this analysis is, how- 
ever, beyond the scope of the present work. 

It should be stressed that, as discussed above, the bal- 
listic RG is redundant for scatterers with l/r potentials, 
like charged impurities or ripples. In this situation, in- 
clusion of interaction does not lead to any essential mod- 
ifications as long as the system is in the ultraballistic or 
ballistic regime. 



VI. SUMMARY 

In this paper, we have analyzed transport properties 
of a graphene sample in the "wide and short" geometry, 
W ^> L, with disorder effects restricted to intra-valley 
scattering. Starting from the clean limit and using the 
transfer-matrix technique, we have analyzed the evolu- 
tion of the transmission distribution P(T) and, in par- 
ticular, of the conductance G and the Fano factor F, 
with increasing system size L. To take the randomness 
into account, we have developed a perturbative treat- 
ment of the transfer-matrix equations supplemented by 
an RG formalism describing the renormalization of disor- 
der couplings. This has allowed us to get complete ana- 
lytical description of the transport properties of graphene 
in the ultraballistic (L <C A) and ballistic (A <C L <C I) 
regimes. We have also constructed "phase diagrams" of 
different transport regimes (ultraballistic, ballistic, diffu- 
sive, critical) for graphene with various types (symme- 
tries) of intra-valley disorder. 
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FIG. 9: Integration contours used to calculate oscillating cor- 
rection to the (a) generating function [Eq. (|A3|) ] and (b) trans- 
mission distribution [Eq. (|A7[) ] . 



apply a saddle-point method. Representing cos(2ueL) as 
a real part of an exponential function and deforming the 
integration contour as shown in Fig. [9H, we get 



8T{z) 



2We 



Re 



i 2 ) 2 du 



2iueL 



V / (1-Z)(1-71 2 ) 3 (1-ZM 2 ) 



(A3) 

The integrand decays exponentially when the contour 
runs far from the real axis. Hence we can estimate the 
value of the integral by expanding the pre-exponential 
factor near the two ends of the contour. We parameter- 
ize these two parts by substituting u = iw and u = 1 + iw 
and obtain the two contributions 

AWe f°° , 3 _ 2weL 3We 



Functional Nanostructures of the Deutsche Forschungs- 
gemeinschaft. The work of I.V.G. was supported by the 
EUROHORCS/ESF EURYI Awards scheme. I.V.G. and 
A.D.M. acknowledge kind hospitality of the Isaac New- 
ton Institute for Mathematical Sciences of the Cambridge 
University during the program "Mathematics and physics 
of Anderson localization: 50 years after" , where a part 
of this work was performed. 



APPENDIX A: OSCILLATIONS AT HIGH 
ENERGIES 

In this Appendix we consider transport properties of 
a clean graphene sample in the limit of high energies, 
eh 1. We will find the next term in the inverse energy 
expansion of the generating function T(z) and the distri- 
bution function P(T), which yields an oscillatory correc- 
tion to the results l|2"Tj) and (p{l"j) . Our starting point is 
the exact expression for the generating function Eq. (|24|) 
that we rewrite using the parameterization 



T(z) = — 

7T 



', du 



cos 2 (ueL) 



sin 2 (ueL) 



(A 1 ) 

Trigonometric functions in the integrand rapidly oscil- 
late. To take advantage of this property, we represent 
the integrand as a sum over Fourier harmonics cos(nueL). 
The first and the second terms of such Fourier expansion 
are 



We 



u du 



1 



iry/T^lJo ^(l-u 2 )^- zu 2 ) 

cos(2ueL) 



2(uyT 



l-u 2 



(A2) 



The first term of the above expression gives the main 
contribution to the generating function, Eq. (|2T|) . The 
next term is suppressed due to oscillations of the inte- 
grand. To estimate this contribution to T{z) we will 



8T h 



IUesin(2eL-7r/4) 
v/2tt(1 - z) 2 



dw 



-2wnL 



We sin(2eL - tt/4) 



8V7?(l-z) 2 (eL)3/ 2 ' (A5) 

We see that the vicinity of u = yields much smaller cor- 
rection, 8T a <§C 8J-b, and hence should be discarded even 
though 8Tb oscillates. The generating function including 
the first oscillating correction is the sum of Eqs. (|27| and 

H), 

K(z)-E(z) 



T(z) = We 



sin(2eL - tt/4) 
zy/T^~z ^ 8^{1- z) 2 (eL) 3 / 2 



(A6) 

The conductance and the Fano factor [Eqs. (|3"3")) and (I34[l ] 
are then calculated by expanding T(z) in small z. 

Let us now calculate the oscillating correction to the 
distribution function (f3"Tj) . The second term in Eq. (|A6|) 
for the generating function does not possess a branch cut 
at z > 1 [this is an artifact of the saddle-point approx- 
imation applied to Eq. (|A3p ], Therefore, we cannot get 
the result by applying Eq. (f3"0"| to the generating function 
[6j) . Instead, we have to use a more general expression 
[2|) . which still possesses a branch cut. Performing the 
analytic continuation of the integrand in Eq. (|A2j) and 
using Eq. ([50]) . we obtain the correction to the distribu- 
tion function in the form 



8P(T) 



2We duu 2 [T - (2 - T)u 2 } cos(2ueL) 

^T 2 Jo ^(\-T)(\-u 2 f{T-u 2 ) 



This integral contains rapidly oscillating function and we 
calculate it using the same method as above. We replace 
cos(2ueL) by exponential, then deform the contour as 
shown in Fig.m>, and estimate the integral in the vicinity 
of u = vT by changing variable u — VT+iw. The result 
of this calculation is 



8P{T) = - 



V2Wesm(2eLVT + tt/A) 

Wesin(2eLVT 



dw 



-2eLh 



IW 

tt/4) 



7r 3/2 T i/4( 1 _ T)VeL 



(A8) 
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Combining Eq. (|A8|) with the main part Eq. (I31j) . we 
have the following distribution function 



K{T)~E(T) sin(2eL\/r + 7r/4) 



P(T) = We 

tt 2 TVT^T 7 r 3 / 2 T 1 /4(i_T)Vei 

(A9) 

The correction to the distribution function, Eq. (| A8[> . 
is not integrable at the point T = 1. This prevents us 
from calculating corrections to conductance and higher 
moments using Eq. (|A9[) . In fact, the result (|A9|) is not 
accurate when T is close to 1. Indeed, expanding the in- 
tegrand in Eq. (|A7|) near u = \/T, we have neglected the 
variation of the factor (1 — it 2 ) -3 / 2 . When T approaches 
1, this neglection is not justified because the singularity 
at it = 1 gets close to the integration contour [see Fig. 
HJd]. The integral in Eq. (|A8|) converges at w ~ 1/eL. 
When typical values of w are of the order of 1 — VT our 
approach fails. Thus we have to impose the condition 



1 



1 

T> — 
eL 



(A10) 



for applicability of the result (|A9|) . This condition also 
ensures that the oscillating correction is smaller than the 
main term in Eq. (|A9j) . 



APPENDIX B: DERIVATION OF THE 
TWO-LOOP BETA FUNCTION FOR RANDOM 
POTENTIAL 

1. The model and universality 

In this Appendix we derive the RG equations for ran- 
dom potential disorder. The beta function is universal 
when it is invariant under small changes in the definition 
of the coupling constants. Generally, the latter depends 
on details of the high-energy part of the spectrum (where 
the dispersion is no longer linear) and hence on the way 
the ultraviolet cutoff of the effective low-energy theory is 
imposed. Therefore, the invariance of the beta function 
with respect to uncertainty in the definition of couplings 
is equivalent to its independence on the RG regulariza- 
tion scheme. 

When only one coupling constant is present in the 
model, the beta function is universal within the two- 
loop accuracy. Indeed, in this case a small change in a 
coupling constant would only change the three-loop and 
higher terms in the beta function. In order to see this, 
one can assume that the beta function for some coupling 
a is known within the three-loop accuracy, 



dc 
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= Aa 2 + Ba 3 + Ca* 7 



(Bl) 



with the coefficients A, B, and C in the first, second, 
and third-loop terms, respectively. Introducing a new 
coupling a' through 



a' = a + Ma 2 +Na 3 , 



and using Eq. (|B1[) . one finds the RG equation for this 
new coupling in the form 

o / 

Aa' 2 +Ba' 3 + (C—AM 2 +AN~BM)a' i . (B3) 
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One sees that the coefficients of the first and second terms 
of the new beta function remain unchanged, whereas the 
coefficient in the last term (third loop) depends on the 
definition of a' . 

When the model contains more than one coupling con- 
stant, already the two-loop RG equations are in general 
non-universal. This can be seen from 



Afa k ai + Bf m a kai a 



kirn , 



<91nA 

a'i = ati + Mfa k a u 
da', 



(B4) 
(B5) 
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him i 



+ (Bf m + 2AfMf n - 2A^M l l n ) a'^a' 



(B6) 

In our model, the RG equation for ao is independent of 
energy e, therefore we can retain the two-loop term. On 
the other hand, the RG equation for energy involves both 
e and ao; hence, we only keep the one-loop term in the 
corresponding scaling function. 

We choose the dimensional regularization (with the 
minimal subtraction) as our RG schemed we consider 
the action in d = 2 — e dimensions (e > 0) and send 
e — > at the end. The model is characterized by two 
constants, a mass parameter m which corresponds to the 
imaginary (Matsubara) energy and the coupling constant 
ao which corresponds to the mean quadratic potential 
disorder strength Eq. (f59"| . 

The renormalized action of the model (known as the 
massive Gross-Neveu model) has the form 



Sr[4>] 



(B7) 

Here m = — it and we have introduced the mass scale /i to 
fix the dimension of coupling ao to zero. The wavcfunc- 
tion i\> is a c?-dimensional spinor in the left/right-moving 
space, cr is a vector of tr Q , which are the generators of 
the (i-dimensional Clifford algebra, obeying 



E 



a a a a = dt. (B8) 



(B2) 



Our goal is to derive the RG equations for ao and m, 
which determine the evolution of these two parameters 
upon increasing the (infrared) scale of the model. This 
derivation closely follows that of Refs. 68.69, where a 
related massless theory was considered. 

We will start with the one-loop calculation, since the 
corresponding integrals and counterterms will be re- 
quired for the two-loop calculation as well. As discussed 
in Sec. [Vl we will discard diagrams with closed fermionic 
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loops, as is appropriate for a system with quenched disor- 
der. (Alternatively, the same result is obtained by using 
the replica trick or the supersymmetry.) 

The quadratic (clean) part of the action yields the bare 
fermion propagator (solid line in diagrams), 



G (0) (p) = (im-ap)- 1 = 



<rp — im 



p z + m 



2 ' 



(B9) 



Dashed lines in the diagrams denote the disorder corre- 
lator 



I\0) = 2ir/i e a . 



(BIO) 



In order to keep track of the two-sided algebra structure, 
we draw the diagrams for vertex corrections with an up- 
per and lower electron line. 

The counterterms are denoted with crossed circles. 
The vertex counterterm is 

ST = 2TT^Sa i <g> 1 (Bll) 

and the self-energy (mass and velocity) counterterm is 

<5E = (-i5m + Sv crp) 1 ® 1. (B12) 

Below we will calculate one- and two-loop diagrams for 
the self energy and the vertex amplitude and construct 
the corresponding counterterms using the minimal sub- 
traction scheme. 

The divergent parts of all the integrals appear with the 
factor c e or c|, where 



and 7 w 0.577 is the Euler-Mascheroni constant. 



(B13) 



2. One-loop RG equations 

We calculate the one-loop integrals with the accuracy 
0(e), because this precision is needed for the two- loop 
calculation later. The following integrals appear in the 
one-loop diagrams of Fig. Q] (we include fi £ to make the 
integrals dimensionless) : 



1 



(27r) 2 - £ p 2 



j2-£ 



p 



PaP/3 



(2tt) 2 - £ (j> 2 



-P 



,2 



0(e) 



- + 0{e) 



(2tt) 2 -^ (p 2 + 



= c e [l + 0(e)]. 



(B14) 
, (B15) 
(B16) 



All other integrals appearing in the one-loop diagrams 
are zero because of isotropy. 

Only two diagrams (Fig. 2^ andSJa) give the first-loop 
corrections to the self energy and vertex, the third and 



fourth diagrams (Fig. [4J: and [4jd) cancel each other [up 
to 0(1)]. More specifically, the diagram (c) on its own, 



(c) = 4 7 r 2 a 2 l c £ 



1 \ - 

i <g> 1 2_^a a ®v 



0(e). (B17) 



would generate a new algebraic structure (corresponding 
to a new disorder — random vector potential). We have 
to calculate it together with its crossed companion, the 
diagram (d): 



(d) = 47T age, 



1 \ - 



+ 0{e). (B18) 



In combination of the two diagrams, the new structure 
is canceled. This happens also in higher loops, where 
all new algebraic structures always cancel. We therefore 
combine diagrams with their crossed versions directly. 
The one-loop corrections read 



E (D 



p=0 



72-e 



P 



1 



(2i:) 2 - e pp 2 + m 2 



-2irim aoc £ 



0(s) 



(B19) 



rW| p=0 = 2x(b) + (c) + (d) 



o 2 2 2e 

87T a Q a a (Tf3 n 



j2-e 



p 



PaPf3 



i a 2 2 e 
l07T a n C £ 



(2tt) 2 ^ {p 2 + m 2 ) 2 
p m 



<l 2 ■• ~ 

J2w)2-e ( p 2 + m 2)2 



O(l) 



+ 0(1) 
(B20) 



The one-loop correction to the self energy is independent 
of external momenta, and therefore, there is no one-loop 
correction to the velocity. Thus no rescaling of the fields 
is required. 

Within the minimal subtraction scheme, the diver- 
gences in Eqs. (|B19|) and (|B2Q|) are canceled by the fol- 
lowing one-loop counterterms: 



o y 'm = . 

e 



5^ 



a 



yielding the one-loop RG equations 



dao 
<91nA 



= 2al 



de 
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2a 2 



ea , 



(B21) 



(B22) 



with A being the real-space running scale (A ~ /i ) and 
energy e = im. 

3. Second loop 

Let us now turn to the second-loop calculation. In 
addition to the integrals appearing in the one-loop di- 
agrams, the further two single-integrals appear in the 
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two-loop calculation: 
m 



(2tt) 2 - £ (p 2 +m 2 ) 3 
d 2_e p m 2 p a p /3 



- 2+ { e) 



(27r) 2 " £ (p 2 + m 



2\3 



(B23) 
(B24) 



The two-loop calculation will only be performed up to 
0(1), since only the divergent parts of the diagrams are 
required for deriving the RG equations. 

The double integrals that appear in the calculation of 
the second-loop diagrams can be reduced to the following 
set: 



,2£ 



d d pd d q 


PaPp 






(27r) 2d {p 2 -( 


- m 2 )(q 2 + m 2 )[(p + q) 2 + m 2 ] 


= <W c\ 




1 




d d pd d q 


Pa qp 






(2n) 2d {p 2 4 


- m 2 )(q 2 + m 2 )[(p - 


- q) 2 + m 2 ] 


= S a p c 2 e 






d d pd d q 


p a qp(p + q)n(p- 





(B25) 



(B26) 



(27r) 2d (p 2 + m 2 )(q 2 + m 2 )[(p + q) 2 + m 2 ] 2 

d „ jd . 
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0(1), 
(B27) 



2s fd d pd d q 



rw-papp 



(27r) 2d (p 2 + m 2 )(q 2 + m 2 )[(p + q) 2 + m 2 } 2 
'1 

e 



lap ' 



•0(1) 



(B28) 



2e fd d pd d q 



m^p a qp 



(2ir) 2d (p 2 + m 2 )(q 2 + m 2 )[(p + q) 2 + m 2 



-- + 0(1) 



(B29) 



The integrals with m 2 and m 4 in the numerator of an 
i ntegr and of the type Eqs. (|B25|) - (|B26j) and (|B27|) - 
(|B29|) . respectively, are not divergent [~ 0(1)]- All other 
integrals appearing in the diagrams can be derived from 
the above set by combining (|B25j) - (|B29|) and/or using 

P *-> q- 



a. Self energy 

Let us calculate corrections to the self energy in two- 
loop order (including diagrams with the one-loop coun- 
terterms). The relevant diagrams are shown in Fig. [10a 
- [TOb . The divergent parts of these diagrams and the 
integrals used in their calculation are given in Table [I] 
The divergent part of the self-energy correction is 



-(2)1 



3i ma 2 . 



lp=0 



(B30) 



TABLE I: Two-loop diagrams. The first column contains 
combinatorial factors and diagram labels according to Fig. 
1101 The integrals relevant for the calculation of a given di- 
agram are listed in the second column. The third column 
contains the divergent part of the diagram. 



diagram 


integrals used 


result 




self energy [in units Air 2 ima^c 2 / e 2 \ 


W 


(|B14()-(|B15(). (|B14()-(|B16() 


-4 + 4e 


(b) 


(|B25(). (|B26() 


-2 


self enerj 


5y with counterterms [in units A-K 2 ima^c e /4tt£ 2 ] 


(c) 


(|B15|) 


4 - 4e 


(d) 


(|B14|) 


8 




velocity 




(e) 


(|B27|). (|B28() 


— 47r 2 crpao c e/£ 




vertex [in units 8n 3 ao)j, £ 


c 2 Je 2 ] 


2x(f) 


(IB15I) Z . (|B15|)-(|B16|) 


8 Sc /Scr\ 
o — oc — \t)C/ 


2x(g) 


(IB27|). (|B29() 


-4 + 4e + (4e) 


4x(h) 


(|B27|). (|B29|) 


8 - 8e - (8e) 


4x(i) 


(1B14()-(|B23(). (|B14|)-(|B24() 


-<8e) 


G) 


(|B15|) Z . (|B15|)-pi6|) 


4 - 4e - (4e> 


4x(k) 


(|B29|) 


-(16e) 


2x(l) 


(|B29() 





2x(m) 


(1B14()-(|B23(). (IB14()-(|B24() 





(n) 


flBlBI)" 


4-2e 


2x(o) 


(|B27l) 


-4 + 4e 


vertex with counterterms [in units 


87r 3 ag/x e c e /47T£ 2 ] 


2x(p) 


(|B15|) 


-16 + 8e + (8e) 


2x(q) 


()B15() 


-16 + 8e + (8e) 


4x(r) 


(|B23(). (|B24|) 


(8e> 


2x(s) 


(|B16|) 


(16e) 


2x(t) 


(|B23|). (|B24|) 






The divergence in E^ 2 ) is compensated by the second- 
order mass counterterm 



3 ma 2 . 



(B31) 



Electron velocity i>o also acquires a correction in two- 
loop order. This correction appears because the second- 
order self energy depends on the external momentum. 
Expanding S^ 2 ^ in small external momentum we obtain 
the diagram Fig. flUe . It equals to 



(e) 



c7£( 2 )(p) 



dp 



erp 4:TT 2 alcl 



p=0 



0(1) 



(B32) 

This divergence is compensated by the velocity countert- 



4e 



(B33) 
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FIG. 10: Two-loop diagrams for the self energy (a-e) and vertex (f-t) corrections. 



b. Vertex 

The two-loop vertex diagrams are shown in Fig. 1 1 OF - 
[TUt . The values of these diagrams are given in the bottom 
part of Table Q] along with the integrals used in their cal- 
culation. According to Refs. [r38l|6 9[ one can disregard the 
mass in the numerator of the electron propagator (|B9j) . 
In order to check this fact, we keep the corresponding 
parts in angular brackets in Table |U Indeed, these con- 
tributions sum up to zero. 

The divergent part of the two-loop vertex correction is 



r (2) 



16d 32c F 2d 



e 2 4tte 2 



^ (__ + _]. (B34) 



This divergence is canceled by the two-loop vertex coun- 
tcrtcrm 



2s 



(B35) 



4. RG equations 

Now we collect the one and two-loop counterterms and 
compose the bare action of the model: 

Sb = I d 2 ~ £ x[^B<TVi>B + rnB$B4>B+na 0B ('ipB'>pB) 2 ] 



I 



= / d'- £ x 



I + I iper\7i/j + m I 1 



an 3ag 



2s 2 



2an 4ag a 2 , \ , 



(B36) 



Parameters of this action are 
^3=^(1 



m b = m 1 



og 
8e 

Oq 

e 



3ag _ ag 
2e 2 4e 



aos = a 0/ u 1 



2ao 

e 



4a 



(B37) 
(B38) 
(B39) 



By construction, the bare couplings tub and aoB do not 
depend on the scale /i, 

pi = 0, |^ = 0. (B40) 
a m /i a In \i 

This determines scaling behavior of renormalized (ob- 
servable) couplings: 



dap 

dm 
d In fi 



-eao — 2a — 2a 



-map 



(B41) 



(B42) 



We express the result in the form of real-space scaling 
with A ~ fi . This amounts to changing the sign of 
the derivatives. Taking the limit e — > and replacing 
mass by the energy (m = — ie), we finally obtain the RG 
equations in the form 

da 



<91nA 

de 
9 In A 



2ag + 2ag, 



= e a + — 



(B43) 
(B44) 



As discussed above [see Eqs. (|B4I) - (|B6|) ]. the second- 
loop term in the RG equation for energy (1B44|) is not 
universal. This is easily demonstrated with the help of a 
small redefinition of olq, 



a' = a + Mag, 



de 
9 In A 



Ma' " 



(B45) 
(B46) 
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On the contrary, the two-loop term in the beta function 
for a remains unchanged as the RG equation (|B43|1 . 
being dimensionless, cannot contain e. In the main text, 
we use only universal one-loop part of the energy beta 
function. 



Using Eqs. (|C6[) and (|C7[) . we obtain the RG equation 
for gi in the form 



dgi 16 a / 1 



(C9) 



APPENDIX C: ANALYSIS OF ONE-LOOP RG 
EQUATIONS 

In this Appendix we will analyze the set of coupled 
RG equations ([75]) . ([77)) and (|75|) assuming that all three 
couplings are nonzero. It is convenient to introduce new 
couplings 



Oil 
Oil 



a±_ + a z , 
a± — 2a z 



(CI) 
(C2) 



In terms of these couplings the one-loop RG equations 
take the form 

dot 2 

g lnA = 2«o + 2a ai + g( 2a i + a 2)(ai - a 2 ), (C3) 
c^cy 2 

— — - = 2a ai + -(ai + 2a 2 )(ai - a 2 ), (C4) 
am A 9 

<9a 2 4 

- = -4a a 2 - -(ai + 2a 2 )(ai - a 2 ). (C5) 

It follows that the couplings ao and ai grow upon the 
renormalization. Furthermore, comparing Eqs. (|C3[) and 
(|C4[) . one sees that the coupling ao increases faster than 
ati. On the other hand, the coupling a 2 decreases, which 
allows us to neglect a 2 in Eqs. (|C3[) and (|C4|) as the first 
approximation: 



2a 



o -r 2a ai 



dap 
dhiA 

9ai 2 2 

_ = 2a a 1 + -a 1 . 



(C6) 
(C7) 



Let us consider the ratio of the couplings, which we de- 
note 



9i 



Oi 

ao 



(C8) 



Assume bare values of the three couplings are of the 
same order, ao ~ a± ^ a z ^ a. Since ao grows faster 
than a\ , for the "first iteration" we retain only ao in Eq. 
(IC91). We find 



1 



1 16 r dA 



x ao(A). (CIO) 

0i (A) 3i 9 J a A 

Further, neglecting ai also in Eq. (|C6|) . dividing this 
equation by ao, and integrating the result over din A, 
we find 



, dA , T . , a (A) 
2/ _ao(A) = ln-^. 
/a A a 



(Cll) 



When a reaches unity (this happens when e < 7), com- 
bining Eqs. (ICIO]) and ([ClT]) . we have 



1 



1 1 



— ■ In — . 

3i(A) ai 31 9 a 



(C12) 



We arrive at the conclusion that at the end of the renor- 
malization the scalar potential becomes the dominating 
disorder: 



Oio 



a± + a z 



a ~l 



9 a 



(C13) 



This result can be refined by retaining a\ in Eqs. 
(|C6|) . (fC7|) . and (|C9)), wh ich yields a correction ~ 
-(3/4) In J In a\ to Eq. (ICl3l . 

Above we assumed that all three couplings are of the 
same order at the ultraviolet scale. The generalization 
onto the case of non-equal couplings is straightforward. 
It turns out that in the ultraballistic regime (when the 
renormalization stops by the largest coupling reaching 
unity) the scalar potential always wins, whatever the ini- 
tial couplings are. 
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